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Abstract 

A new model for stock price fluctuations is proposed, based upon an 
analogy with the motion of tracers in Gaussian random fields, as used 
in turbulent dispersion models and in studies of transport in dynamically 
disordered media. Analytical and numerical results for this model in a spe- 
cial limiting case of a single-scale field show characteristics similar to those 
found in empirical studies of stock market data. Specifically, short-term 
returns have a non-Gaussian distribution, with super-diffusive volatility, 
and a fast-decaying correlation function. The correlation function of the 
absolute value of returns decays as a power-law, and the returns distribu- 
tion converges towards Gaussian over long times. Some important char- 
acteristics of empirical data are not, however, reproduced by the model, 
notably the scaling of tails of the cumulative distribution function of re- 
turns. 

PACS numbers: 05.60.-k, 89.65.Gh, 05.40.-a, 02.50-r 

1 Introduction 

Random fluctuations in stock market prices have long fascinated investors and 
mathematical modelers alike. Although the investors' hopes of accurately pre- 
dicting tomorrow's share price appear to be in vain, models which limit them- 
selves to statistical characteristics such as distributions and correlations have 
had some success. Bachelier's classical model treats the stock price S{t) 
as a random walk, leading to the conclusion that the distribution of prices is 
Gaussian. Samuelson 2, instead describes the log-price 

x{t) = Hs{t)/sm (1) 

as a random walk, and therefore concludes that the stock price S{t) should 
have a log-normal distribution. This model remains in common usage, despite 
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the shortcomings listed below, not least because it permits the derivation of 
an equation for pricing options and other financial derivatives, e.g., the famous 
Black-Scholes equation ,3j. Nevertheless, empirical evidence from stock mar- 
ket data indicates that the random walk model inadequately describes many 
important features of the stock price process. The following stylized facts are 
accepted as established by these studies [31 0] : 

(i) Short-term returns are non-Gaussian, with 'fat tails' and high central 
peaks. The center of the returns distribution is well fitted by Levy distribu- 
tions [Ej. Recent studies indicate that the tails of the returns distribution 
decay as power-laws 6 . As the lag time increases, the returns distribu- 
tions slowly converge towards Gaussian; recent analysis of the Standard & 
Poor 500 index (S&P500) estimates this convergence is seen only on lags 
longer than 4 days 0. 

(ii) The correlation function of returns decays exponentially over a short 
timescale, consistent with market efficiency. However, the correlation 
function of the absolute value of the returns shows a much slower (power- 
law) decay 

(iii) The volatility (standard deviation of returns) grows like that of a diffusion 
(random walk) process, i.e. as the square root of the lag time, for lag 
times longer than about 10 minutes. However, higher frequency (shorter 
lag) returns demonstrate a super-diffusive volatility, which can be fitted 
to power-laws with exponents found to range between 0.67 and 0.77 [Hl[7]- 

(iv) The distribution of stock price returns exhibits a simple scaling: in js] a 
power-law scaling of the peak of the returns distribution P(0) with lag 
time is shown to hold across many magnitudes of lag times. The exponent 
of the power-law is approximately —0.7. 

Most of the references cited here examine data from the Standard & Poor 500 
index (S&P500), but other international markets are found to behave similarly 



The predictions of the random-walk model can be investigated using the 
random differential equation 



which yields the log price x{t) from each realization of the random function 
Uuj{t). The classical random walk follows from taking u^it) to be a Gaussian 
white-noise process with zero mean and (auto)correlation function [5]: 



where 6{t) is the Dirac delta function, and the angle brackets denote averaging 
over time or over an ensemble of realizations. White noise is the formal deriva- 
tive of a Wiener process, which connects Q to rigorous methods for stochastic 




(2) 



a^S{t - t') 



(3) 
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differential equations Pj- Note that throughout this paper all stochastic pro- 
cesses are assumed (unless explicitly stated otherwise) to have zero mean, zero 
drift, and to be statistically stationary. 

Equation (|2Jl is easily solved, and leads to results which do not agree with 
the empirical facts (i)-(iv) listed above. Defining the return on the stock price 
S{t) at time t as the forward change in the logarithm of S{t) over the lag time 
A: 

r-A(t) = x{t + A)-x{t) 

= \n S{t + A) ~ In Sit), (4) 

it is found that the white-noise case leads to probability distribution functions 
(PDFs) for TA which are Gaussian at all lags A, in contrast to (i). The lack of 
memory effects in white-noise leads to correlation functions which decay imme- 
diately to zero, unlike (ii). The volatility grows diffusively for all lag values, and 
so has no super-diffusive range, while the scaling of P{0) follows a power-law 
different from that found for actual prices, see (iv) above. 

A natural generalization of the white-noise model again takes a form similar 
to @ 

(5) 

but with Uc{t) now being a non- white or colored process, incorporating some 
memory effects. For the particularly simple case of Gaussian noise, it can be 
described fully by its correlation function 

{u,it)uM = a^R{t - t'). (6) 

Here R{t) decays from a maximum of unity aX t — Q to zero as t — > oo, and 
is the variance of the process Uc- Gaussian noise is attractive because it 
appears naturally as the limiting sum of many independent noise sources under 
the central limit theorem [H|, and because of its analytical tractability. Indeed 
equation lO is again easily solvable for Gaussian Uc, but leads once more to 
Gaussian-distributed returns at all lag times, in contrast to (i). Moreover, the 
correlation of absolute returns does not decay more slowly that the returns 
correlation, as required by (ii). 

These failings of simple differential equation models have inspired many at- 
tempts to reproduce the empirical facts (i)-(iv) using various stochastic models. 
As a sample of just a few of these, we list the truncated Levy flights model 
UnmH, arch and GARCH models [12 US], non-Gaussian Ornstein-Uhlcnbeck 
processes ^3], and a model based on a continuous superposition of jump pro- 
cesses ^Id- The first three of these were reviewed in jE]- An ideal model 
should reproduce all the experimental facts (i)-(iv), by giving a simple picture 
of the stochastic process underlying the stock price time series, but no model 
has yet been found which fulfills all these criteria |17) . 

In this paper we examine the results of generalizing the random walk models 
(jSl and (0) to motion in a random field: 

'^^u{x,t). (7) 
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Here u(x,t) is a Gaussian random field with zero mean, which is fully described 
by its correlation function 

{u(x,t)u{x',t')) = a^Q{x - x' ,t~t') (8) 
= a^S{x-x')R{t-t'). (9) 

For simplicity we take the field to be homogeneous and stationary; the factor- 
ization in (PI of the correlation into separate time- and ^-correlations is for ease 
of exposition only, and more complicated inter-dependent correlations can be 
considered. The averaging procedure in equation (jS)) is over an ensemble of 
realizations of the random field u, or using a uniform measure over all possible 
{x, t) values. 

Equation implies that the log-price x{t) changes in a random fashion, but 
the rate of change is randomly dependent on both time and the current value 
of X. Note that random-walk models lO and © may be recovered by taking 
^(a:) = 1 and choosing the time correlation R{t) appropriately. We will show 
that the x-dependence in the noise term yields qualitatively different results to 
the random walk models, and that many of the empirical observations (i)-(iv) 
of stock markets may be reproduced by Q) . 

The x-dependent noise term in Q is motivated by recent studies in turbu- 
lent dispersion ^21 and transport in dynamically disordered media j^. These 
investigations employ equation (|7|), but with a rather different interpretation: 
x(t) is the position vector (in 2 or 3 dimensions) of a passive tracer, and u{x,t) 
is the velocity vector field which transports the tracer. The tracer is called pas- 
sive because it is assumed not to affect the velocity field by its presence, and 
so u{x^ t) is a prescribed random field. An instructive example is the motion of 
a small buoy on the surface of the ocean 120]: the two-dimensional vector x{t) 
then describes the position of the buoy (by its longitude and latitude, say) while 
the tracer is moved by the ocean waves according to Q , with the Gaussian field 
u{x, t) providing an approximate description of the ocean wave field. 

The distribution of such tracers resulting from motion in a Gaussian field is 
known to be non-Gaussian, at least over short to intermediate timescales ^1 
I21| . Crucial to understanding this effect is the difference between the Eulerian 
velocity u{x, t), i.e., the velocity measured at the fixed location x at time t, and 
the Lagrangian velocity 

v{t) = uix{t),t), (10) 

which is the time series of velocity measurements made by the tracer itself as 
it moves through the random field l22j . When there is no x-dependence in the 
velocity field u, the Eulerian and Lagrangian velocities coincide, but otherwise 
a Gaussian Eulerian velocity field can yield non-Gaussian Lagrangian processes. 
The relationship between Eulerian and Lagrangian random processes is an active 
research area, especially in the case of a compressible (non-zero divergence) 
velocity field which is of interest here [231 1241 1^ . 

Our model of the log- price uses a scalar "position" x{t) and "velocity" u{x,t) 
rather than the vector-valued analogues in the turbulent dispersion problems 
described above. Nevertheless, the concept of Eulerian and Lagrangian velocities 
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carries over to the one-dimensional case, and so we will refer to x{t) as the "price 
tracer", with its "velocity" u{x,t) related by equation Q. The fact that the 
Eulerian velocity field u{x,t) is Gaussian allows us to obtain analytical results, 
and may be justified as a consequence of the central limit theorem applied to the 
sum of many independent forces on stock prices. In this paper we will not further 
pursue a microeconomic justification [201 El for Q, but rather concentrate on 
demonstrating that it may yield some of the observed statistical properties of 
empirical stock market data. 

The remainder of this paper is structured as follows. In section 2 the conse- 
quences of the x-dependence of the noise in equation O are examined through 
some simple examples, and a special limiting case (single-scale field) is high- 
lighted. The statistics of the Lagrangian velocity H10|l in the single-scale case 
are found analytically in sections 3 and 4, and are related to the statistics of the 
returns ta- Numerical simulations of equation (0 are considered in section 5, 
and a simple method for generating Monte-Carlo time series in the single-scale 
case is found. Results from such simulations are presented in section 6, and 
compared to the stylized facts (i)-(iv), as presented in |S] and [S]. Section 7 
comprises of discussion of the results and directions for future work. 



2 X-dependent velocity fields 

The concept of x-dependence in the velocity of the log- price x{t) was introduce 
in section 1. To gain some intuitive understanding of the effect of such x- 
dependence, we first return to the random walk model with colored noise Uc{t), 
independent of x, as in ((SJ and (jSJ: 

dx 

— = uJt), with {uJt)uJt')) = a^Rit - t'). (11) 
at 

For clarity, we consider in this section a specific form for the time correlation 
R(t): 

R{t) ^cxp(-^y (12) 

where a characteristic decay time T has been introduced. Figure 1 shows the 
contours of a single realization of the x- independent Uc{t), negative values of the 
velocity being denoted by shaded areas and dotted contour lines. Time increases 
along the horizontal axis, with the stock log price x measured on the vertical 
axis. The motion of a stock tracer with initial condition a;(0) = is shown also; 
note the negative rate of change of stock price in the shaded regions (negative 
velocity), and the positive rate of change in the regions where u > 0. Because 
there is no x-dependence in 1)11(1 . the velocity field depends only on time, and 
so the contours are vertical lines. By contrast, a simple x-dependent velocity 
field of the type proposed in lO has a non-trivial x-correlation function S{x) in 
equation Q; for example in Fig. 2 we use 



^(^)=exp(-^), (13) 
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where I defines a correlation scale {I ^ oo recovers the x-independent case) . The 
contours in Fig. 2 display a more complicated structure than in Fig. 1, and this 
is reflected in the observed motion of the stock prices. The three parameters a, 
T and I characterizing the random field may be combined into one dimensionless 
quantity 

and if the time, log price and velocity are re-scaled using these parameters: 

t = Ti 
X — Ix 

u = au, (15) 
then the equation of motion (0 becomes 

e^=u{S:,i). (16) 

Note that the re-scaled velocity field u has unit variance. In Fig 2. the stock 
price motions resulting from solving (|16|l are shown for e values ofe=l,e = 0.1, 
and e — 0.01. When e is large, the length scale / of the randomness in x is much 
larger than the effect aT of the random variation in time — this is analogous 
to the case of weak space dependence in turbulent dispersion problems studied 
in using perturbation theory. In the case of e ^ 1, small deviations of 
the stock price distribution from a Gaussian shape are found, with kurtosis less 
than 3. However, the studies of empirical data mentioned in section 1 indicate 
that the kurtosis of stock returns is much larger than 3, and so the weakly 
space-dependent theory does not seem relevant here. 

The opposite limit, i.e., e ^ 1 appears more promising. In Fig. 2, we observe 
that the stock prices for parameter values e = 0.1 and e = 0.01 display some 
interesting generic features. Constrained by the equation of motion to move 
downwards in shaded regions (where the velocity u is negative), and upwards 
in unshaded regions, the stock tracers tend to become trapped near curves of 
u = 0, with u positive below them and negative above them. These stable 
positions disappear when the zero- velocity manifold 'folds over', and the stock 
tracer is then driven by equation Hlfi|l to move quickly to another slow manifold. 
A complete understanding of the dynamics of the system as e ^ requires a 
singular perturbation approach |2H| to equation II16I) , but the important feature 
for our work is the following of the zero-velocity manifold by the stock tracer, 
punctuated by fast jumps. 

An analytical approach to the general e ^ case has not yet been found, 
but some interesting features arc highlighted by choosing a particularly simple 
form for the a;-correlation function: 

S{x) = cos(fcx), (17) 

where fc is a constant. This correlation function corresponds to a so-called 
'single-scale' Gaussian velocity field, which (as shown in [^) can always be 
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written in the form 

u{x, t) = f{t) cos{kx) + g{t) sin(fcx), (18) 

where f{t) and g{t) are independent Gaussian random functions of time, each 
with zero mean, variance a^, and a given correlation function R{t): 

{f{t)f{t')) = {g{t)g{t'))^a^R{t-t'). (19) 

An example of motion in such a single-scale field is shown in Fig. 3. Note the 
periodicity in the vertical direction, and the lack of jumps. In the remainder of 
this paper we will concentrate on this simplified example and show that many 
of the statistical characteristics of the resulting stock returns are remarkably 
similar to those observed in empirical data. 

3 Lagrangian velocity and PDF 

In this section we consider the statistical characterization of motion along the 
It = curves of the single-scale random field (|18|l . according to the equation of 
motion 0, and with the Gaussian random functions f{t) and g{t) defined as 
in (|19|) . Again, the analogue with tracer motion on the surface of a turbulence 
ocean proves instructive. The Eulerian velocity at a point x at time t is defined 
by the random field H18() . Each tracer particle moving through this field will, 
however, experience its own time history of velocity variations; this is called the 
Lagrangian velocity of the tracer v{t). The crucial point in our use of random- 
field models of stock price motion is that the Eulerian field may, as in H18|l , have 
Gaussian statistics, while the Lagrangian velocity (that is, the velocity 'felt by 
the tracer') need not be Gaussian. Some intuitive notion of why this might be 
may be gained from considering the fact that the velocity at space-time points 
(x, t) picked at random from the Eulerian field has a Gaussian distribution; 
on the other hand, the dynamics of the tracer particles mean that they are 
more likely to cluster in regions of low velocity, with high-velocity transitions 
between, meaning both zero values of v and high values are more likely than for 
a Gaussian distribution. This argument holds for any random field, but as we 
show in this section, a quantitative description may be derived in the case of a 
single-scale field (|17|l for the e ^ limit. 

In the limit of small e, we saw in the previous section that the stock tracer 
motion is confined to curves with u = Q. While the Eulerian velocity on these 
curves is zero by definition, the Lagrangian velocity is non-zero, as the tracer 
moves to remain on the u = curve. In fact, the Lagrangian velocity felt by 
the particle trapped on such a curve is given by 



du 




dt \ 




du 




dx 


\u^0 



(20) 



with the derivatives evaluated on the curve u = 0. For the single scale field 
1)18(1 . an explicit expression for the Lagrangian velocity may be found in terms 
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of the functions /(t) and g{t) and their derivatives: 



vit) = '-i:^. (21) 



When the Lagrangian velocity v{t) is known, the equation of motion for the 
tracer is simply 

dx , ^ , , 

(22) 

with solution ^ 

x{t) = x{Q) + I v{t')dt'. (23) 

Similarly, the return on the stock over a time lag of A is given by 

r^{t) = x{t + A)~x{t) 

v{t')dt'. (24) 

Clearly the statistics of the returns follow from the statistics of the Lagrangian 
velocity, indeed for short time lags it might be expected from (I24II that 

TAit) « vit)A, (25) 

so that the distribution of returns, for instance, can be related directly to the 
distribution of the Lagrangian velocity. We shall see later that H25|) is not fully 
accurate for large values of u, but the Lagrangian velocity distribution will still 
prove extremely useful. 

To find the probability distribution function P{v) of the Lagrangian velocity, 
we first consider the cumulative probability that the Lagrangian velocity v at 
a given time is greater than some chosen value V. According to 1)21(1 . v > V if 
and only if 

/' > ^ kV +^ = F, (26) 
9 9 

where we use the symbol F to denote the right-hand-side. Since /, g, f and g' 
are independent Gaussian random variables at any single moment in time, we 
denote their PDFs by P/, Pg etc., and the cumulative probability of the event 
(EHl is 

dg df dg' / df'PgPfPg,Pf,. (27) 

-oo J —OO J —OO JF 

The PDF of v then fohows from 



V=v 



/oo pca POO Qp 

dg df dg'PgPfPg, PiAf,^^ — 
-OO ^ — oo <J —OO ^* 



(28) 



V=v 
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and the fact that 

Notmg that the variances of /' and g' both equal ~a^R"{0) — /t^, where 
we introduce the notation tq for the timescale defined by the initial radius of 
curvature of the function R: 

To = , ^ (30) 

the integrations over the f-g plane may be performed by using polar coordinates 
/ = rcosO, g — rsin6', yielding 



/oo /'OO 1 

dg / d/P,P/^=a=/FT7exp 



2a2 



(31) 



2[1 + k^T^V^]^ 

The PDF (|31|l of the Lagrangian velocity P{v) is symmetric in v, and so has 
mean zero. Note that the tails of P(w) decay as \v\~^ for large \v\, so the variance 
of V does not exist. Also, (pTT|l effectively contains only one free parameter fcro, 
and is otherwise independent of the choice of correlation function R{t) for / and 

9- 



4 Correlation of Lagrangian velocity and returns 

The (normalized) correlation of the stock returns with lag A, at time t, is defined 
by averaging the product of the returns at times separated by r: 

{rA{t)rAit + T)l ^^^^ 



{rAm ' 

and for stationary returns becomes independent of t . The correlation is nor- 
malized by the variance (rA(t)^) of the returns at lag A, which is known as 
the squared volatility In this section we derive analytical results to show 
how the volatility and the returns correlation for a single-scale field with e — > 
depend on the time correlation function R{t) of the Eulerian field, and on a and 
k. 

Because of the dependence (|24|l of the returns on the Lagrangian velocity, 
the correlation of the returns can be written in terms of the correlation function 
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of v: 



rt+A pt+T+A 

{rA{t)rAit + T)) = I dh dt2 {v{h)v{t2)) 

It + T 

(•A pT+A 

dh J dt2 {V{ti)v{t2)) 

cA pT+A 



J dh J dt2L{h-t2) (33) 

/ iA~h)[Lih-T)+L{h+T)]dh 
Jo 



(34) 

where we have used the stationarity of the returns, and defined the Lagrangian 
velocity correlation function 

L{h ~ h) = {v{hHt2)) . (35) 

Note the Lagrangian velocity is assumed to be a stationary process — this has 
not been proven to follow from Eulerian stationarity in the general case |23) . 
but appears to hold in numerical simulations. We will proceed to calculate L{t), 
and then use (|34|l to find the volatility and returns correlation. We note here 
that when the separation time r is much larger than the lag A, the returns 
covariance ()34|l is approximated by 

{rAit)rA{t + T)) ^ A^Lir), (36) 

consistent with (|25|l . 

The correlation function of the stationary random process v is written as an 
eight-dimensional integral, using the definition (|21() of the Lagrangian velocity: 

L{t) = {v{0)v{t)) 

1 /" / /l.9l - g'lfl /2g2 - .92/2 p ,r H f H'.n I I l\ 

= 12 / / f2 , 2 „2 , 2 Pi2(Jl,h,hJ2)Pl2(gi,9l:92:g2)- 

1^ J J Ji + 5i J2 + 92 

(37) 

Here the eight-fold integrals are over /, g, /' and g' evaluated at each of the two 
times t = (with subscript 1) and t = t (subscript 2). The random functions 
/ and g are independent of each other, but the joint probability functions for / 
and its derivative at different times is found from the Gaussian joint PDF 



J'i2(w)--i-^(DetA)^exp 



1 T» 

-— w Aw 
2 



(38) 



for the vector of arguments, w — (/i, /{, /2, /2), with a similar expression for g. 
The matrix A is defined by 

A^' = i^^wj) (39) 
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and when written out in full: 

1 R{t) -R'{t) 

^ ^ R'{t) 1 ^ 

-R"{t) -i?"(0) 

Using H38|) and H40() in (|37|l enables us to calculate the Lagrangian correla- 
tion L{t). We briefly describe here the steps in the calculation of the multi- 
dimensional integrals. First, the four integrals over the variables /{, /2, g'2 
are performed - the dependence of the integrand on each of these variables is 
linear, so the Gaussian integrals may be calculated straightforwardly. For the 
remaining four variables, we transform to polar coordinates in the (/i,5i) and 
if 2, 92) planes: 

/i=rcos0, gi—rsinO, /2 = /0COS(/>, g2—psa\4>. (41) 

The 4> integral is then trivial, and by integrating over 9, r, and p, the final result 
emerges: 

m = 2^5^ {Rm"it) ~ R'm log [1 - Ritf] . (42) 

Equation 1)42(1 gives an explicit formula for the Lagrangian velocity corre- 
lation in terms of the time correlation R{t) of the single-scale Eulerian field. 
When the function R{t) is given, L{t) follows immediately from (|42l) . and hence 
the volatility and returns correlation can be calculated using (|34|l . Before ex- 
amining the results for the returns, we consider the limiting forms of L{t) for 
small and large arguments. 

Assuming that R{t) is sufficiently smooth near t = to allow the definition 
of To as in equation , a small-time expansion of yields 

^^^^^"2fc^^°s(f ) 

Note this diverges as t approaches zero, so the variance of the Lagrangian ve- 
locity does not exist, in accord with our conclusion at the end of section 3. 

For large time arguments, we are particularly interested in power-law decays 
of correlations, and so consider how the decay of R(t) influences that of L{t). 
Taking a power-law decay with exponent — /3: 

R(t) ^ff^ as t ^ 00, (44) 

the corresponding Lagrangian correlation is found from H42|) to decay as 



Lit) - -Jkt'^^'^ ast^oo. (45) 



The negative sign here indicates that the Lagrangian correlation approaches 
zero from below at large times, being negative whenever R{t) has a power-law 
decay. 
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The squared volatility is found by setting t = in H34|) : 



vol^(A) = (rAitf) = 2 / (A - t)L{t)dt. (46) 
Jo 

For short lags, with A ^ tq, it follows from 143() that the volatility may be 
approximated by 



1 A 



, A2 

3 - log ( 

'0 



A' 



T, 



O — (47) 



Note that while this expression limits to zero as A vanishes, it does not follow 
a simple power law in A because of the logarithmic term. Although empirical 
stock volatilities have been fitted with power-laws |S1 17] , we show in the next 
section that (|47|l can also match the data quite well. 

For large time lags A, empirical data shows a linear growth in vol^ with A. 
From H46|l . we can find an expression for the rate of change of vol^ with A: 

d 



dA 



vol^(A) = 2 /" L{t)dt, (48) 
Jo 



and hence conclude that the squared volatility grows linearly with A for large 
lags if the integral 

poo 

L{t)dt (49) 

has a finite positive value. Assuming a power-law form of R{t) for large argu- 
ments, we use the result H45II for the asymptotic form of L{t). It immediately 
follows that the integral H49|l is finite if (3 > —1/2, and hence at large lags the 
squared volatility grows linearly with A in this case. 

In empirical stock data, correlations of nonlinear functions of the returns of- 
ten display significantly slower decay than the correlations of the returns them- 
selves. One quantity of interest is the correlation of absolute returns, which, 
using H25|l . may be related to the correlation of absolute Lagrangian velocities 
(for short lags) as 

(\rA{t)\ \rAit + r)|) - {\rAitf\) « A^ {(HtHt + r)|) - (|«(t)^|)} . (50) 

Following similar procedures as in the calculation of L{t), it can be shown that 
the large-t asymptotics corresponding to a power-law R{t) ~ give a corre- 
lation of absolute Lagrangian velocities scaling as 

{\vio)vm-{y\)-t-"'- (51) 

In contrast to H45|l . this correlation remains positive at large times; note too 
that its rate of decay is slower than that of (|45|l . 

Before proceeding to numerical simulations, we briefly summarize the theo- 
retical results. We have obtained the PDF of the Lagrangian velocity in (|31|l . 
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and related the Lagrangian velocity correlation L{t) to the given Eulerian corre- 
lation R{t) in H42I) . This permits estimates of the volatility as a function of lag. 
Finally, we expect a power-law scaling for the absolute returns correlation when 
R{t) has a power-law decay. In the following sections we perform numerical 
simulations to confirm and extend these results. 

5 Numerical simulations 

To extend the analytical results found in the previous section, we consider the 
implementation of numerical simulations of motion in random fields. Random 
fields u{x, t) may be generated using standard techniques [3T], with the ordinary 
differential equation j?)) being solved using Runge-Kutta methods. Averages 
may be calculated over an ensemble of realizations, or over a long time series 
x{t) in a single realization, to closely mimic the statistical analysis of S&P 500 
data performed in E] . 

A random field with zero mean and correlation function Q may be generated 
using a superposition of random Fourier modes |31|: 



with the amplitudes An and _B„ chosen from independent Gaussian distributions 
of zero mean and variance . The w„ and fc„ are chosen from distributions of 
random numbers chosen so as yield the correlation (PJ. Specifically, the the 
LOn are chosen from a distribution shaped as the Fourier transform of R(t), 
with the distribution of the fc„ being the Fourier transform of S{x). Thus, 
Gaussian distributions with zero mean and unit variance are used for the tj„ 
and fc„ to generate the random field of Fig. 1, as the Fourier transforms yield 
R{t) = cxp(— and S{x) — exp(— a;^/2) as required in equations (|12|l and 



For the special case of a single-scale field considered here, the fc„ are all i/c, 
and the field may be written in the simpler form (|18|) . see [30]. Accordingly 
we only require a method for generating the random functions of time f{t) and 
g{t), and this is easily derived from ifS^ : 



with a similar formula for g{t). The An, Bn are chosen as above. We are 
especially interested in the effects of power-law correlations R{t), so we choose 
the ujn in H53fl from the Gamma distribution |31| : 




(52) 




(53) 



G{u) 



'UJ' 



(54) 
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Here T{/3) denotes the usual Gamma function. The Fourier transform of G{uj) 
is the correlation function R{t) of /: 



y2 



cos 



/3tan ^ 



T 



(55) 



which decays as R{t) ^ t^^^ for t > T. 

The methods described above are sufficient to simulate motion in a random 
field as described by or equivalently by the rescaled equation H16|l . and 
indeed this is how the stock log curves in Figs. 1 to 3 are generated. However, 
in the idealized limit e ^ considered in the previous sections, the tracer closely 
follows the contour u = 0. In this limit it is therefore not necessary to explicitly 
solve the differential equation to determine x{t) in a single-scale field, as this can 
be determined from the implicit equation u{x, t) = 0. From (|18|l we immediately 
obtain 



. ^ 1 -1 

x{t) — — tan 
k 



m 



(56) 



for motion along the u ~Q contour. Given the expression H53|l for f{t) and g{t), 
equation H56|l then generates the time series x(t). Some care must be taken to 
ensure continuity of x{t) near times when g[t) = 0, but overall this method is a 
very efficient way to create a model time series for x{t). 



6 Numerical results 

In this section we report the results of employing the algorithm (|56|l to generate 
a sample time series x{t) corresponding to a single-scale field, with Eulerian time 
correlation function given by H55(l . Having decided upon a single-scale field and 
the e ^ limit, so that the stock tracer follows contours of u{x,t) = 0, our 
remaining choices to fit empirical stock data is rather limited. We are free to 
choose the time correlation function R{t) and the spatial scaling factor fc; note 
the standard deviation a of the random field has disappeared because we have 
take the limit e 0. 

Our choice for R{t) is motivated by the empirical results such as those re- 
ported in |6J, and by our analysis in previous sections. We have seen in section 
4 that the Lagrangian velocity correlation becomes negative, with a power-law 
decay exponent of —2/3 — 2 when R{t) has a power-law decay exponent of —f3. 
On the other hand, the correlation of absolute values of the velocity remains 
positive for large t, with a decay exponent —2/3. Studies of the corresponding 
returns correlations in |S] have concluded that the correlation decays exponen- 
tially quickly to the noise level, whereas the absolute returns correlation exhibits 
a slow decay with power-law exponent of approximately —1/3. Matching this 
to our result (|51|l motivates us to choose a power-law form for R{t) as given in 
l|55(l . taking the value oi (3 — 1/6. The time-scaling T is chosen to be 10 min- 
utes, to approximately match the volatility results of [S]- Having thus defined 
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the function R{t), we are left with only the single parameter k. This is chosen 
to be 500, to approximately fit the volatility results (see below) to those in 

Each realization consists of 1.6 x 10^ values of x{t), modelling the log price 
at 1-minute intervals, to mimic the S&P500 data set used in [S]. Averaging 
is over time within each realization; we show results from different realizations 
only when statistical scatter is evident in the single-realization results. From 
the series for x{t), it is straightforward to calculate the series of returns over 
integer time lags A 

rA{t)=x{t + A)-x{t). (57) 

and to calculate various statistical properties of the returns. 

The volatility 

vol(A) = VirAm (58) 

found by numerical simulation is plotted as a function of the lag A in Fig. 4 
(filled circles). The parameter k is chosen as fc = 500, in order to closely match 
the empirical A = 1 volatility value as shown in Fig. 3(c) of (Hj. Also shown in 
Fig. 4 is the analytical result H4t)|l . plotted as a dotted line, which confirms that 
the numerical method closely reproduces the exact results of section 4. The 
dashed line in Fig. 4 is the small- lag approximation 147|) . and shows that the 
super-diffusive region of volatility (for lags under 10 minutes) is well- matched by 
the logarithmic-corrected quadratic power law in A. Analysis of the empirical 
data frequently leads to fitting the super-diffusive volatility with a power-law, 
with exponents in the range 0.67 6 to 0.77 |7j. At longer lags (A > 10), random 
walk behavior (vol ~ a/A) is observed in all studies, and reproduced in our 
model as shown following equation H49(l . The solid lines in Fig. 4 show power- 
law scalings, with exponents 0.77 (for A < 10), and 0.5 (for A ^ 10), which are 
found in [7| to match the S&P500 volatility. Clearly these scalings also match 
the results of our model quite well, indicating that the model predictions behave 
similarly to the empirical volatility curve. Indeed, our result 147() predicts that 
a log-corrected quadratic fit should be superior to the power-law fits used in 
[HI [7]. Testing this prediction requires higher frequency analysis of the empirical 
market data. 

The probability distribution function (PDF) of the one- minute returns, ri{t), 
is plotted in Fig. 5. Open circles are from numerical simulations, and are well- 
fitted by the Lagrangian velocity PDF H31|l . up to high returns values. Note the 
returns are normalized by their standard deviation a = vol(l) = 2 x 10^^. The 
Gaussian distribution with this standard deviation is plotted with the dotted 
line: note the higher-than-Gaussian central peak of the model results, and the 
fatter tails. Fig. 5 is remarkably similar to Fig. 2 of which shows the 
one-minute returns PDF from the S&P500. Indeed, in |S] it is demonstrated 
that the empirical results are well-fitted near the center of the distribution by 
a Levy stable distribution with index a = 1.4, and so in Fig. 5 we also show 
this distribution (dashed line), with scale factor chosen to match the peak value 
of the numerical distribution. Comparison with Fig. 2 of 5 shows that our 
numerical simulation returns have a very similar distribution to the S&P500 
data, at least within ±10cr. 
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In Fig. 6 we consider the correlation function of the one-minute returns over 
a separation time r, normahzed by the volatihty: 



{n{t)n{t + T)) 



(59) 



As this quantity has both positive and (smaU) negative values, we plot its ab- 
solute value on a log-linear scale. Statistical scatter causes some uncertainty 
in these numerical results, so the results of two different realizations (circles 
and triangles, respectively) are shown. According to H36|l and (|45l) . the returns 
correlation is negative for r 3> 10 minutes, which is confirmed by plotting the 
analytical result (|34|1 : this is shown as a solid line where the correlation is pos- 
itive, and a dashed line where it is negative. The numerical simulation results 
match the analytical form well for r < 10 minutes, but degrade in quality when 
the correlation is negative. It is possible that more advanced numerical simu- 
lation methods such as Fourier- Wavelet could alleviate this problem. Note 
filled symbols show positive numerical correlation values, with open symbols 
for negative values. The rapid decay of the correlation from r = to r = 10 
minutes (see dotted line showing exp(— t/3) for comparison) is reminiscent of 
that found in S&P500 data, see Fig. 3(a) of and Fig. 8(a) of |2H|- Unlike 
our model results, however, the empirical correlation does not exhibit negative 
values but instead reaches a "noise level" around 3 x 10"'^, where it flattens out 
without discernable structure. Our model suggests that it might be fruitful to 
search for evidence of negative correlations in the empirical data, which may 
currently be screened by the noise level. 

The correlation function of the absolute value of one-minute returns (nor- 
malized to unity at zero separation) is 



From empirical data, this is known to have a slow power-law decay with r, with 
exponent of approximately —1/3, see Fig. 3(b) of 6 and Fig. 8(b) of [2H|- The 
analytical result (|51|l for the present model indicates that a similar power-law 
form holds, indeed this is what motivated our choice of parameter /3 = 1/6 in 
H55|) . Fig. 7 shows numerical results from two separate realizations (circles and 
triangles), with the solid line showing a power-law decay with exponent —1/3. 
The numerical correlations are all positive, but show a slight deviation from 
the expected power-law scaling. This may be related to the poor representation 
of the negative values of the returns correlation in Fig. 6. Comparison with 
empirical data (Fig. 3(b) of [Hj and Fig. 8 (b) of 28 ) also shows that the 
model's correlation decays faster than that of the S&P500 data for r < 10 
minutes. 

The non-Gaussian distribution of returns over longer time lags is examined 
in Figs. 8 and 9, using the methods applied to S&P500 data in Figs. 6 to 8 of 



{\n{t)n{t + r)\)-{\n{t)\f 
ilnm) - {\riit)\f 



(60) 
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The normalized returns are defined by dividing by the volatihty: 



9A{t) 



rA{t) 
vol(A) ' ■ 



(61) 



and the moments 



(l5A(i)|™) 



(62) 



of the numerical simulations are plotted in Fig. 8(a) as a function of m, for lags 
of 1, 16, 32, 64, and 128 minutes. The extremely non-Gaussian distribution at 
A = 1 minute relaxes to a shape which remains stable up to A « 128 minutes — 
note the moments for lags from 16 to 128 minutes are almost indistinguishable 
in Fig. 8(a). This conclusion is supported by the cumulative distribution of the 
normalized returns shown in Fig. 8(b). This is very similar to the behavior ob- 
served in empirical data, see Fig. 6 of [Hj , although the model does not produce 
power-law scaling of the distribution tails (except for the region with exponent 
-2 for lag A = 1 minute). The returns distributions eventually converge to 
a Gaussian distribution, see the corresponding plots for lags of 390, 780, and 
1560 minutes (1 to 4 trading days) in Fig. 9. This convergence to Gaussian is 
somewhat faster that that observed in the S&P500 in |6), where non-Gaussian 
moments and cumulative distributions persist until approximately 4 days. Our 
corresponding estimate for the model persistence time (based on the moments 
in Fig. 8(a)) is 128 minutes, an order of magnitude smaller that the S&P500 
value. 

The scaling of the peak P(0) ofthe PDF of stock price returns S{t+A) — S{t) 
was examined in |H] for S&P500 data. A power-law scaling of P{0) with lag was 
discovered, with exponent —0.7. In Fig. 10 we mimic this study, generating the 
time series S{t) from x{t) in order to calculate the stock price return. Results 
are shown from two separate realizations, along with a line indicating Gaussian 
scaling of exponent —1/2. It is clear that, unlike the S&P500 returns, the 
model i-'(O) does not scale as a single power-law with A, although it appears to 
approach the Gaussian scaling at the highest lags. 

7 Conclusion 

We have proposed a new model for the fluctuations of stock prices, in which 
the rate of change of the log price is randomly dependent upon both time and 
the current price, see equation (jj)). This is a very general modelling concept, 
with the attractive feature that Gaussian Eulerian fields may give non-Gaussian 
Lagrangian statistics in the observable data, i.e. the time series of returns. As 
noted in section 5, Monte-Carlo simulations of equation 10 can be used to find 
the model predictions for any random field. In this paper we focus on a special 
case, the e — > limit of a single-scale random field, for which exact analytical 
results are obtainable, and for which numerical simulations are particularly 
efficient. 

Our main analytical results are the expression 142(1 for the Lagrangian veloc- 
ity correlation L{t) in terms of the given Eulerian field correlation R{t), and the 
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quadrature formulas (|34|l and H46|) giving the correlation function and volatility 
of the returns. Numerical simulations of stock price time series are efficiently 
performed using equation H56|) . and allow us to examine features not amenable 
to exact analysis. 

We find that several of the important empirical stylized facts (i)-(iv) listed 
in the Introduction are reproduced by our model, using the correlation function 
H55(l for computational convenience. The parameter values k and T are chosen 
by comparing the model's volatility to the data in Fig. 3(c) of [HI. With reference 
to the stylized facts listed in section 1, the results of the single-scale random 
field model may be summarized as follows. 

(i) Model returns are non-Gaussian, with fat tails and the center of the dis- 
tribution well-fitted by the Levy distribution used in [S]. The tails of 
the cumulative distributions of model returns do not, however, have the 
power-law scaling found in jH]. The normalized model returns distribu- 
tions exhibit a slow return toward Gaussian, on timescales an order of 
magnitude longer than the characteristic time of R{t). However, these 
timescales are an order of magnitude smaller those found for the conver- 
gence to Gaussian in the S&P500 returns in 0. 

(ii) The correlation function of model returns decays quickly (similar to expo- 
nentially) over a short timescale. Following the fast decay, the correlation 
becomes negative, with magnitude decaying more slowly to zero. Although 
the fast decay for lags under 10 minutes is very similar that found for the 
S&P500 in 1^, the empirical correlation function does not reach negative 
values, rather remaining at a constant "noise level" . Whether this noise 
level could be masking negative correlation values as predicted by the 
model is a question requiring further extensive analysis of market data. 

(iii) The volatility of the model returns grows diffusively for lag times longer 
than 10 minutes. The super-diffusive volatility at higher frequencies (shorter 
lags) is actually due to a logarithmic correction to power-law growth (|47() , 
but for the range of lags examined in C] it matches well to a super- 
diffusive power-law. Higher frequency data is needed to determine which 
form best fits the actual stock market volatility. 

(iv) Our model does not have a simple scaling of P{0) with lag, in contrast to 
the case of S&P500 prices as demonstrated in Fig. 1 of fSJ. 

In summary, the single-scale random field model yields time series which have 
many (but not all) of the important properties of empirical data. The model 
gives an appealingly intuitive picture of the causes of non-Gaussian behavior in 
markets, and provides a simple algorithm for generating time series with many 
market-like features. We anticipate that this algorithm will prove particularly 
useful in Monte-Carlo simulations to determine prices of options and derivatives 
when non-Gaussian returns are particularly important |21 El ■ 

The concept of modelling stock prices by motion in a random field is quite 
general, and suggests many directions for further work. Among these are: 
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(a) Examining the consequences of such models for the pricing of options and 
derivatives, and comparison with the standard Black-Scholes (random- 
walk-based) results. 

(b) Further market data analysis, to determine if there is any evidence of the 
negative correlation or log-corrected volatility scaling discussed in (ii) and 
(iii) above. 

(c) Finding a microeconomic basis for the random field model Q, possibly 
based on the actions of multiple independent agents |26[ I77| . 

(d) Numerical simulations of iQ for more general field correlations S{x) and 
different e values, to determine which of the results detailed here are uni- 
versal features of motion in random fields, and which are specific to the 
single-scale, e ^ case analyzed here. 
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t/T 

Figure 1: Motion of a stock tracer in a random field which depends only on 
time, as in ((SJ. Shaded areas and dotted contour lines indicate negative values 

of Uc{t). 




Figure 2: Motion of stock tracers in a random field which depends on both time 
and X, as in (0. The correlation functions are (|12|l and Ijl^ll . and the values of 
e are: e = 1.0 (dash-dotted), e = 0.1 (dashed), and e = 0.01 (solid). 
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Figure 3: Motion of stock tracers in a single-scale random field. See Fig. [5] for 
description of line types. 




Figure 4: Volatility of returns as a function of lag time A. Symbols show the 
results of numerical simulation; dotted line is the analytical expression (I46|) : 
dashed line is the small- A approximation H47|l : solid lines show power-laws with 
exponent 0.5 (large lags), and 0.77 (small lags), as fitted to empirical data in 

El- 
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Figure 5: Probability distribution function of 1-minute returns from numerical 
simulation (symbols). Dashed line is the Levy distribution fitted to empirical 
data in 0; solid line is the Lagrangian velocity PDF JSIJ; dashed line is a 
Gaussian distribution with same variance as the data. 
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Figure 6: Magnitude of the correlation function of 1-minute returns, as a func- 
tion of separation time r. Numerical results from two realizations are shown 
(circles and triangles), with the analytical form H34|l (solid and dashed lines). 
Negative correlation values are shown with open symbols and dashed line. The 
dotted line shows exp(— r/3) for comparison. 
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Figure 7: Correlation function of absolute 1-minute returns, as a function of 
separation time r. Numerical results from two realizations are shown (circles 
and triangles). The solid line is a power-law scaling as H51|l. 



25 




Figure 8: (a) Moments of the normalized retmns, for lag times A = 1, 16, 32, 
64, and 128 minutes, (b) Corresponding cumulative distribution functions, with 
the dashed line showing a power-law with exponent -2. 
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Figure 9: (a) Moments of the normalized returns, for lag times A = 390, 
and 1560 minutes, (b) Corresponding cumulative distribution functions. 
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Figure 10: Scaling of the peak P(0) of the stock returns PDF with lag time A. 
Solid line shows a Gaussian scaling, with exponent —0.5. 
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